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ABSTRACT 

We take a sample of 24 elliptical and lenticular galaxies previously analysed by the 
SAURON project using three-integral dynamical models created with Schwarzschild's 
method, and re-analyse them using the made-to-measure (M2M) method of dynamical 
modelling. We obtain good agreement between the two methods in determining the 
dynamical mass-to-light (M/L) ratios for the galaxies with over 80% of ratios differing 
by < 10% and over 95% differing by < 20%. We show that (M/L) M 2M ~ (M/L) Sc h- For 
the global velocity dispersion anisotropy parameter 5, we find similar values but with 
fewer of the made-to-measure models tangentially anisotropic by comparison with their 
SAURON Schwarzschild counterparts. Our investigation is the largest comparative 
application of the made-to-measure method to date. 

Key words: galaxies: elliptical and lenticular - galaxies: kinematics and dynamics - 
galaxies: structure - methods: N-body simulations - methods: numerical 



1 INTRODUCTION 

Within the field of galactic and stellar dynamics, it has 
become common practice to model kinematic observa- 
tions of a galaxy in order to interpret the observations 
1 and to understand better the underlying dynamical struc- 
tures withi n the galaxy. Within this modelling arena, the 
method of ISchwarzschiidl <jl979h has been heavily devel- 
oped and deployed with over 500 citations from other 
pape r s (for example, iRix et al.l Il997l , Ivan den Bosch et al.l 
20081 . IJalali fc Tremaind l201lh By comparis on, the made- 
to-me asure method (M2M) formulated bv lSver &; Tremainel 
(1996) is less well-known but is no less capable, and 
has be en the subject of growing interest recently (for ex- 
ample, Jourdeuil & Emsellem 2007, de Lorenzi et al. 2007, 



de Lorenzi et al.1 120081 . iDehnenl 120091 . lLong fe Maol I201C . 
Das et al.l 120 111 ). Both methods achieve their objectives by 



weighting a system of particles / orbits and superimposing 
them to reproduce the galactic observations. The key differ- 
ence is that in Schwarzschild's method a library of orbits is 
first created and then weighted, whereas in the M2M method 
the orbit weights are determined dynamically as the parti- 
cles are being orbited . Other methods exist w hich, while not 
directly derived from lSver &; Tremainel {1996), seek to tailor 



the kinematics of a system of particles to match the kine- 
matics of a galaxy, for example Rodi onov et al. (2009). 

In this pape r we compare the mad e-to-measure method, 
as descr ibed by | Long &; Maol (1 20 10), and Schwarzschild's 
method. Cappellari et al. (2006) use Schwarzschild's method 
to determine the mass-to-light ratios for a selection of galax- 
ies observed with the SAURON0 integral- field spectrograph. 
These same galaxies are re-analysed using the M2M method 
and the resulting ratios compared. The galaxies comprise a 
mixture of elliptical and lenticular galaxies covering both 
fast and slow rotators, and including both edge-on galaxies 
and galaxies inclined to the line of sight. As an extension to 
the mass-to- light e xercise, we calculate t he global anisotropy 
parameters as in I Cappellari et all (|2QQ7h and again comp are 
the results. Earlier pap ers (for example, iDas et al.ll201ll or 
Ide Lorenzi et al.1 12009) have used the M2M method effec- 
tively with individual galaxies. To our knowledge, this pa- 
per is the first to use the M2M method with a larger sample 
of galaxies and is the first to compare directly the results 
achieved with those from using Schwarzschild's method. 

In section [2] we describe the M2M method, and in sec- 
tion [3] its application to the SAURON galaxies. Sections 2] 
and [5] cover respectively the mass-to-light determinations 
and the global anisotropy parameters. We draw the activi- 
ties to a conclusion in section [6] As might be expected, we 
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refer heavily to the published SAURON material for data 
values. We do not however cover in detail any theory from 
the SAURON material unless there is some specific point 
to be made in relation to the M2M method. Unless other- 
wise stated we adopt the same modelling assumptions as 
ICappellari et al.1 <|2006h . 



2 THE M2M METHOD 
2.1 Outline 

In brief, the M2M method is concerned with modelling stel- 
lar systems and individual galaxies as a system of test par- 
ticles orbiting in a gravitational potential. Weights are as- 
sociated with the particles and are evolved over many or- 
bital periods such that, by using these weights, observa- 
tional measurements of a real galaxy are reproduced. We 
expect that the weights themselves will have converged in- 
dividually to some constant value. It is natural to relate the 
particle weights to the luminosity of a galaxy and then to 
consider how the galaxy's surface brightness and luminosity 
weighted kinematics could be generated using the particle 
system. 

In the next section, based on lLong &; Mad (|2010l ). we 
set out the theory underlying the M2M method. 



2.2 Theory 

For a system of N particles, orbiting in a gravitational po- 
tential, with weights Wi, the key equation which leads to the 
weight evolution equation is 



J k 



(i) 



where % 2 5 S and d are all functions of the particle weights 
w = (wi, — -,wn)] t \s time; and fi and e are positive param- 
eters. The equations governing weight evolution over time 
come from maximising F(w) with respect to the particle 
weights (dF/dwi = Vz) and rearranging terms to give 
equations of the form 



-€WiG(l 



(2) 



dt ' 

The overall rate of weight evolution is controlled by e. The 
precise form of the function G(w) depends on the constraints 
d and is illustrated later (equation [9j). The process be- 
ing applied to % 2 is one of regularised, parameterised con- 
strained extremisation. 

The x 2 term in F arises from assuming that the proba- 
bility of the model reproducing a single observation can be 
represented by a Gaussian distribution and then construct- 
ing a log likelihood function covering all observations. For 
K multiple observables, we take % 2 m the form 



(3) 



where A& are small, positive parameters whose role is ex- 
plained in section 13.71 



and 



Vk,j(w) - Y k j 



(4) 



(5) 



where Ykj is the measured value of observable k at position 
j with error (7kj, and ykj (w) is the model equivalent of Ykj. 



j(w) = y^WiK k j(ri,Vi)6(i e kj) 



(6) 



where Kkj(ri,Vi) is the kernel for observable k evaluated 
at position j for a particle with position Vi and velocity Vi. 
5(i G k,j) is a selection function and signifies that only par- 
ticles which contribute to observable k at position j should 
be included in the calculation of ykj. We have listed the 
kernels required for this paper in section 
The entropy function S in F is 



S(») = -5>ln(£: 



(7) 



where mi is taken as the initial value of a particle weight (in 
practice, we take nra — 1/N). S is used for regularisation / 
smoothing purposes with the amount of regularisation being 
controlled by the parameter ji. The derivative term dS/dt 
indicates that over time we require the particle weights, and 
thus iff, to be constant. As demonstrated in lLong &; Mad 
(|20ld ). the term behaves as the constraint dS/dt — 0. 

The functions C% in F are additional constraints to be 
included in the maximisation of F. In this paper, we use 
only one such constraint which is that we require the model 
luminosity to match the luminosity (L) of the galaxy being 
modelled, that is Lwj = L, or more concisely — 1. 

We therefore take 



Ci 



As 



E 1 



(8) 



where A SU m is a positive parameter. 

Given the definitions of % 2 5 S and d and noting that 
for the purposes of this paper we do not use regularisation 
(/i = 0), G(w) from equation [2] can now be written 



G(W) Afc ^ ,J ' A fc)J - + Asum I ^ Wj-1 



(9) 



where Kkj(ri,Vi) has been abbreviated to Kkji and we 
have assumed that, for all observables, 1 particle contributes 
only at 1 position j. 

Finally, model observables (and thus particle weights) 
are subject to noise as the numbers of particles contributing 
to the observables vary. This noise is suppressed by replacing 
Akj in G(w) by an exponentially smoothed version Akj 
given by 

j t K k j = a(A k j -A k j) (10) 

where a is a small positive parameter. The smoothed Akj 
can be used to calculated a smoothed version ykj of the 
model observable, 
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Ykj + (Jk,j^k,j 



(ii) 



and 



3 SAURON M2M MODELS 
3.1 Galaxies and observables 



The g alaxies which we model are as in ICappellari et all 
(|2006b but with NGC 221 omitted since it is not part of 
the SAURON data release. The galaxies are listed in Table 
[1] together with the properties which are relevant to M2M 
modelling. As indicated earlier, the galaxies comprise a mix- 
ture of elliptical and lenticular galaxies covering both fast 
and slow rotators, and including both edge-on galaxies and 
galaxies inclined to the line of sight. The galaxies also ex- 
hibit various core features, for exa mple kinematically dis - 
tinct cores or counter rotating cores (Emsellem et al. 2004h. 

The inclinations a nd distances (di s tance modulus) to 
the galaxies are as per ICappellari eta D (|2006T ) Table 1. We 
have not attempted to use M2M modelling to determine 
the inclinations. Within a M2M model, we employ Carte- 
sian axes such that the positive x-axis points towards the 
observer and the y — z plane represents the galaxy's on 
sky projection. We align the galaxies' photometric major 
axes to the mode l y- ax is utilising position angles taken from 
ICappellari et al.1 (|2007T ). 

We take kinematic data for the galaxie s from the 
SAURON data release (|Emsellem et al.M2004l ). The data 
available are the line-of-sight mean velocity, velocity disper- 
sion and the /13 and /14 Gauss-Hermite coefficients, all taken 
from a truncated Gauss-H ermite expansion of the line -of- 
sight velocity distribution rtvan der Marel fe Franx|[l993h . 



losvd(v) 



exp(-^g orm /2) 



l + ^h n H n { 

n—3 



Vn 



(12) 



where H n is the Hermite polynomial of degree n and the 
normalised velocity v norm is defined as 

v — v 

^norm — (1^) 



where v and a are the line-of-sight mean velocity and veloc- 
ity dispersion respectively. 

The /15 and h§ Gauss-Hermite coefficients used in 
ICappellari et al.1 J2OO6) are not available in the data release. 
We assume therefore that h$ — h§ — with a measurement 
error of ±0.3 (M. Capellari private communication). We do 
not model a galaxy's mean line-of-sight velocity and veloc- 
ity disper s ion d irectly but i nstead model hi = h<i — as in 
iRix et al.1 (|l997l ). Following livlagorrian fe Binnevl (| 19941 ). we 
calculate the measurement errors Ah\ and A/12 as 

Av 

Ah! 4- (14) 



V2a 



and 



Ah 2 = ^ (VI2/14 - y/2) 



2a 



(15) 



where Av and Act respectively are the measurement errors 
in the mean line-of-sight velocity v and velocity dispersion 
a. If we require the model mean line-of-sight velocity v m or 
the model line-of-sight velocity dispersion cr m , we calculate 
them as 



2ah 2 



12/i4, m — \[2 



(17) 



(16) 



where the hj, m are the expon entially smoothed model hi 
values. Ide Lorenzi et al.l (|2009h use a similar approach in 
their M2M models of NGC 3379. It is possible to calculate 
and the Gauss-Hermite coefficients by fitting Gauss- 
Hermite series directly to the end of modelling run particle 
data, provided sufficient particles are available to populate 
the velocity histograms necessary to the fitting process. The 
approach above, using smoothed model values, avoids the 
need to run the M2M models with large numbers of particles. 

We put the SAURON kinematic data through a clean- 
ing process (section 13. 3[) before subtracting the systemic 
galactic velocity from the mean line-of-sight velocity, sym- 
metrizing the data and converting it to units appropriate to 
our M2M modelling (distances in effective r adii, time in 10 7 
years) . We take the systemic velocities from lEmsellem et al.l 
(2004) and assume they are subject to a measurement error 
of 10%. The usual error propagation rules are applied. 

The observables in our M2M models are thus 

(i) surface brightness, 

(ii) Gauss-Hermite coefficients hi to he- 
Values for the kinematic observables are as described above, 
and surface brightness is calculated from the mult i- Gaussian 
expansions of the galaxy's surface brightness (see section 
13. 4[) . For modelling purposes, we assume a 10% relative error 
in surface brightness values. Unless explicitly stated, lumi- 
nosity density is not used in our M2M models to constrain 
the luminous matt er distribution ( s ee sec tion 13. 9[) . 

Similarly to ICappellari et al.1 ([2006), we perturb the 
line-of-sight particle coordinates by a 'point spread function' 
before binning any model data to cr eate the model observ - 
ables. We use the seeing values from lEmsellem et al.l |2004) 
Table 3 and implement the function as a circular Gaussian 
distribution. 



3.2 Voronoi Tessellation 

To achieve a pre-determined signal to noise, the SAURON 
observatio ns were adaptively b i nned and p rocessed, as de- 
scribe d in ICappellari &; Copinl (|2003h and lEmsellem et al.l 
(2004), resulting in a centroidal Voronoi tessellation. The 
kinematic data in the SAURON data release are presented 
in the context of that tessellation. The M2M models use the 
same Voronoi tessellation for determining model kinematic 
observables. The Voronoi cells are the bins used for accumu- 
lating particle kinematic data as part of the construction of 
the model observables. 

Within our M2M implementation, the Voronoi bins are 
represented by a 2D 'tree' with a surrounding convex hull. 
Particles are binned by first determining whether they are 
inside the convex hull and, if they are, performing a 'nearest 
neighbour search' of the 2D tree to identify the bin required. 
The areas of the bins (needed for model observable calcula- 
tions - see section 13. 6ft are calculated using a Monte Carlo 
approach. As an example, Figure [1] shows the positions of the 
Voronoi centroids (the data measurement points) for NGC 
3156 together with their convex hull. The SAURON data re- 
lease contains no information on the extent of the outermost 
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Table 1. Galaxy sample and measured parameters 



Galaxy 


Type 


Fast 


Distance 


Re 




i 


PA 


Seeing 


V S yst 


(M/L) Schw 


(M/L) M 2M 






rotator 


Mpc 


arcsec 


km s 


deg 


deg 


arcsec 


km s 


I-band 


I-band 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 


NGC 524 


S0+(s) 


VPS 


23.34 


51 


235 


19 


48.4 


1.4 


2353 


4.99 


6.39 


NGC 821 


E6? 


VPS 


23.44 


39 


189 


90 


32.2 


1.7 


1722 


3.08 


3.37 


NGC 2974 


E4 


VPS 


20.89 


24 


233 


57 


43.5 


1.4 


1886 


4.52 


4.60 


NGC 3156 


SO : 


VPS 


21.78 


25 


65 


67 


49.4 


1.6 


1541 


1.58 


1.46 


NGC 3377 


E5 - 6 


VPS 


10.91 


38 


138 


90 


41.3 


2.1 


690 


2.22 


2.22 


NGC 3379 


El 


yps 


10.28 


42 


201 


90 


67.9 


1.8 


916 


3.36 


3.67 


NGC 3414 


SOpec 


no 


24.55 


33 


205 


90 


179.9 


1.4 


1472 


4.26 


4.56 


NGC 3608 


E2 


no 


22.28 


41 


178 


90 


79.3 


1.5 


1228 


3.71 


3.73 


NGC 4150 


S0°(r)? 


yps 


13.37 


15 


77 


52 


147.0 


2.1 


219 


1.30 


1.26 


NGC 4278 


El - 2 


yps 


15.63 


32 


231 


45 


16.7 


1.9 


631 


5.24 


5.61 


NGC 4374 


El 


no 


17.87 


71 


278 


90 


128.2 


2.2 


1023 


4.36 


4.65 


NGC 4458 


EO - 1 


no 


16.75 


27 


85 


90 


4.5 


1.6 


683 


2.28 


2.32 


NGC 4459 


S0+(r) 


yps 


15.70 


38 


168 


47 


102.7 


1.5 


1200 


2.51 


2.76 


NGC 4473 


E5 


yes 


15.28 


27 


192 


73 


93.7 


1.9 


2249 


2.91 


3.12 


NGC 4486 


EO - 1+pec 


no 


15.63 


105 


298 


90 


158.2 


1.0 


1274 


6.10 


7.05 


NGC 4526 


SAB0°(s) 


yes 


16.44 


40 


222 


79 


112.8 


2.8 


626 


3.35 


3.26 


NGC 4550 


SB0° : sp 


yes 


15.42 


14 


110 


84 


178.3 


2.1 


413 


2.62 


2.78 


NGC 4552 


EO - 1 


no 


14.93 


32 


252 


90 


125.3 


1.9 


351 


4.74 


5.01 


NGC 4621 


E5 


yes 


17.78 


46 


211 


90 


163.3 


1.6 


456 


3.03 


3.07 


NGC 4660 


E 


yes 


12.47 


11 


185 


70 


96.8 


1.6 


1089 


3.63 


3.85 


NGC 5813 


El - 2 


no 


31.33 


52 


230 


90 


134.5 


1.7 


1947 


4.81 


4.69 


NGC 5845 


E : 


yes 


25.24 


4.6 


239 


90 


143.2 


1.5 


1474 


3.72 


4.34 


NGC 5846 


EO - 1 


no 


24.21 


81 


238 


90 


75.2 


1.4 


1710 


5.30 


5.38 


NGC 7457 


SO-(rs)? 


yes 


12.88 


65 


78 


64 


125.5 


1.3 


845 


1.78 


1.52 



Colum n (1): NGC number. Column (2): Hubble type from Emsel lem et al.l (120041). Colum n (3): Fast/slow rotator from Cappellari et al.l 
(l200d). Column (4): Distanc e calculated from distance modulus in lcappellari et al.l (120061 ) Column (5): Effective (half-light) radius 
from lCappellari et al.l (|2006|). Column ( 6): Velocity dispersion within effec tive radius from Cappellar i et al.l (|2006h . Column (7): 
Inclination from Cappellari et all |2006). Column (8): Positi on angle from [Cappellari et al .1 (120071 ). Column (9): Seeing from 
lEmsellem _el_alJ 120041) Colu mn (10): Systemic velocity from Emselle m et al.l (|2004t ). Column (11): Schwarzschild mass-to-light ratio 
from I Cappellari et al. (2006). Column (12): M2M mass-to-light ratio - section 0] 



Voronoi bins, and, as a consequence, the model observable 
calculations for such bins are biased inwards. This is not 
considered a significant issue as the number of such bins is 
small. 

Voronoi bins are not used with surface brightness. In- 
stead we employ a 16 x 16 pola r grid with radial d ivisions 
pseudo-logarithmic as described iLong fe Maol (|2010h . 

To give an indication of the number of observations be- 
ing modelled, of the 24 galaxies, NGC 5845 has the least 
number of Voronoi bins with 164 and NGC 4486 the greatest 
with 2112 bins. The total number of observations is therefore 
1240 for NGC 5845 and 12, 928 for NGC 4486. 



3.3 Data cleaning 

We have subjected the SAURON data release to a number 
of checks to eliminate suspect kinematic data 

(i) 'nan' (not a number) value check 

(ii) zero value check, for example error fields should not 
be zero 

(iii) positive value check, for example error fields should 
not be negative 

(iv) small value check (< 10 -5 ) 

(v) record sequence check - a record is nagged if a data 



NGC3156 




Figure 1. Measurement points (X$, Y$) are in blue and the con- 
vex hull is indicated by the red dashed lines. The bounding rect- 
angle is used in calculating the bin areas. The circle has the maxi- 
mum possible radius for a circle to lie within the convex hull, and 
is used for performance reasons as part of the algorithm to deter- 
mine whether or not a particle is inside the hull or not. Distances 
are in units of effective radii R e for the galaxy. 



© 2012 RAS, MNRAS OOO.ITJfHl 



Made-to-measure galaxy models II 5 



field in one record has the same value as in the previous 
record 

Note that not all checks are applied to all data fields. 

In total, 13 galaxies were found to have suspect data as a 
result of the exercise, and consequently we have not used the 
associated Voronoi data bins in the modelling process. The 
bins are not actually removed from the Voronoi tessellations 
but marked as 'not in use'. 

As identified in Emsellem et al. | (|2004h . the NGC 5846 
data contain contamination from a foreground star and a 
companion (NGC 5846A). This north and south contamina- 
tion has, as far as possible, been removed by simply deleting 
all the Voronoi bins where the magnitude of the y-coordinate 
position is greater than 3.5 kpc. 

3.4 Gravitational Potentials 

All the galaxies are modelled as axisymmetric galaxies 
with their gravitational potentials calculated from depro- 
jecting the mult i- Gaussian exp ansions (MGEs) of the ir sur- 
face brightness record ed in ICappellari et all ([2006T ) , and 
iKrainovic et al.1 (|2005h for NGC 29 74. The mult i- Gaussian 
expansion technique is described in Emsellem et al.l (1994) 
and is not repeated here. The galactic kinematic and photo- 
metric symmetry axes are assumed to align. No M2M mod- 
ellin g of the galaxi e s with the axes not aligned, as discussed 
in Emsellem et al. I (|2QQ7h for example, is undertaken. 

In our M2M implementation, to avoid multiple numer- 
ical integrations, the MGE potential and associated accel- 
erations are pre-calculated and held on 1500 x 1500 R — z 
interpolation grids with bilinear interpolation used between 
grid points. We employ OpenMP0to accelerate production 
of the grids. 

We augment the MGE potential with a central black 
hole modelled as a Keplerian potential. The mass of the 
black hole Mb h is calculated using t he Mbh - <J relationship 
as described in iGiiltekin et HI (1 20091) (a is the bulge velocity 
dispersion) . For consistency with lCappellari et alj(|2006h we 
do not include a dark matter component in the potential. 

Orbit integration is performed using the standard 2 nd 
order interleaved leap frog method with an adaptive time 
step. The duration of our M2M models is inversely propor- 
tional to the dynamical time of the galaxy being modelled 
and is numerically, approximately 400 divided by the dy- 
namical time with a minimum of 200 units. We use the 
formula for dynamical time in iBinnev fe Tremaind (2008) 
§2.2.2 and calculate it at the half light radius of the model. 
The size of a galaxy model is 3 times the maximum disper- 
sion in the galaxy's surface brightness mult i- Gaussian ex- 
pansion. In practice, the model sizes range from 5 to 18 
effective radii depending on the galaxy. 

3.5 Particle initial conditions 

In setting the initial spatial and velocity coordinates for par- 
ticles, two issues need to be addressed. The first is how to 
handle global rotation of the galaxy, and the second, how 
to handle core features. In both cases, we may choose to 

2 http: / / openmp.org 



take no explicit action and allow the M2M method to at- 
tempt to weight the particles such that the observables are 
reproduced. Alternatively, we may use our knowledge of the 
galaxy's features and set the initial particle conditions ac- 
cordingly. The first approach, taking no explicit action, is 
inefficient in the use of particles (consider the case of many 
particles orbiting in the opposite sense to any global rotation 
- the method will lower the particles' weights to reduce the 
particles' influence on the model observables). We therefore 
discount the first approach and adopt the second. 

We employ two schemes for setting the initial condi- 
tions. In the first, we set the initial spatial positions of the 
particles to approximate the luminosity distribution gener- 
ated by deprojecting the galaxy's surface brightness MGE. 
Creation of the velocity coordinates follows a 3 stage pro- 
cess, 

(i) use the velocity dispersions created by solving the 
semi-isotropic Jeans equat ions for an axisymmetric system 
(|Binnev fe Tremaind l2008h to provide initial values for the 
velocity coordinates, 

(ii) set the global rotation sense for a prespecified frac- 
tion of the particles to align with the rotation sense of the 
SAURON observations, and 

(iii) for the particles inside the SAURON measurements 
convex hull, adjust the coordinates to approximate the mea- 
sured line-of-sight velocity by sampling from a Gaussian 
distribution formed from the measured line-of-sight veloc- 
ity and dispersion. 

The effect of stage (iii) is to reproduce (approximately) 
in the particle initial conditions any core features in the 
SAURON velocity measurements. Determination of the 'pre- 
specified fraction' in stage (ii) is not yet an automated pro- 
cess. The fraction is determined iteratively by comparing 
visually the SAURON velocity contours with the particle 
equivalents and then adjusting the fraction as necessary. As 
an example, we show the resulting velocity contours for NGC 
2974 in Figure [2] 

For the second scheme, a grid-less energy, angular mo- 
mentum and pseudo th ird integral system similar to that in 
ICappellari et al.l (|2006h is adopted to determine the initial 
conditions for the particle orbits. We then modify the initial 
conditions as in stages (ii) and (iii) above. Unless otherwise 
stated, by default, all modelling runs are performed using 
this scheme. 

We determine the number of particles to be used in 
our M2M models by examining the particle distribution 
across the Voronoi bins. We adjust the number of particles 
such that the minimum mean number of particles per bin 
is greater than 60. We find that between 10 5 and 5 x 10 5 
particles are required depending on the galaxy. Ideally the 
minimum number of particles per bin should be used but 
this would require much larger numbers of particles (> 10 6 ) 
and correspondingly more computing resources. Note that 
at any one time during a modelling run significant numbers 
of particles (> 30%) are outside of the convex hull and are 
not contributing to reproducing the kinematic observables. 

3.6 Kernels 

The kernels are similar to those in lLong fe Mad (120101 ) and 
we list them below. 
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Figure 2. For NGC 2974, the left panel shows the symmetrised measured mean line-of-sight velocity and the right panel, the line-of-sight 
velocity from the particle initial conditions (at the start of modelling). The velocity units are km s — 1 . 



K 



(i) luminosity density 

L 

(ii) surface brightness 



(18) 



(19) 



(iii) mean luminosity-weighted Gauss-Hermite coefficient 



h n 



(20) 



where L is the luminosity of the galaxy being modelled, Aj 
is the area of the bin at position j and Ij the target surface 
brightness, Vj is the bin volume, v\\^ is the line-of-sight ve- 
locity for particle z, and H n is the Hermite polynomial of 
degree n. The normalised velocity v n0 rm, ji is denned as 

_ V \\,i - Vj / 01 x 
^norm, ji — \^^) 



where Vj and aj are the measured line-of-sight mean veloc- 
ity and velocity dispersion respecti vely. We normalise the 
Herm ite polynomials as described in Ivan der Marel &; Franxl 
(|l993h appendix A. 



3.7 Parameter setting 

In this section we describe how we set the values of the 
various parameters within the M2M method. As indicated 
earlier, we do not use regularisation for any of the modelling 
runs and so we take \i — 0. For exponential smoothing we 
take a common approach across all the galaxies and set a — 
0.05. The e parameter controls the overall rate of weight 
evolution and we set e = 10 -4 initially. We may alter it 
later if we find that we are not achieving a xl P er degree of 
freedom value of O(l) for a modelling run. Finally, we set 
Asum — 10 4 . 

The role of the observable Xk parameters (see equations 
[3] and [9]) is to help balance the weight evolution equation 
across all the observables. The equation contains terms of 
the form 



a k,j 



(22) 



The Akj component reflects how well the model is repro- 
ducing the measured observations and is not examined fur- 
ther. The Kk,ji/(Tk,j component varies, by several orders of 
magnitude, between observables and between positions j for 
a single observable. By running a M2M model for a short 
period of time (5 dynamical time units), we are able to un- 
derstand how the Kk,ji/crk,j values are varying. We take the 
modal value of Kkji/crkj (found by binning logarithmically) 
and set Xk such that 



A/c ^-|modal — 10. 

&k,j 



(23) 



Similarly to e, we may adjust the value of Xk if we find that 
we are not achieving a xl P er degree of freedom value of 
O(l) for a modelling run. 

For the 24 galaxies we analysed, involving some 168 
Afc's, approximately 25 % of the Afc's required adjustment. 
ICappellari et all |2006) noted that reproducing the Gauss- 
Hermite coefficient Ha proved problematic. Based on the Xh 4 
values we achieve, we do not have an equivalent issue with 

/l4- 



3.8 Computer performance 

Our M2M software has been parallelised using the Message 
Passing Interface (MPI) with the parallelisation being based 
around a star netw ork with a s i ngle c entral controlling node. 
We reported in lLong &; Mad ([20101 ) that the impl ementa- 
tion was high ly scaleable and oth ers, for example iDehnenl 
(2009|) and|de Lorenzi et al. (2007), have reported similarly. 
This position remains true for a low number of observables 
and measurement points. However given the number of mea- 
surements available for the galaxies we have analysed (see 
section |3J"[) , we find that the scaleability is reduced. This 
reduction in our case is due to the overheads of handling 
data packet fragmentation particularly on the central node 
of the network. Increasing the packet size will recover some 
of the reduction. For larger M2M models, it may be appro- 
priate to introduce a layer of nodes whose primary role is to 



© 2012 RAS, MNRAS OOO.ITlfHl 



Made-to-measure galaxy models II 7 



act as data concentrators. We have yet to investigate either 
of these schemes. 



3.9 Miscellany 

As in dicated earlier, this paper builds on lLong &; Mad 



((20101) (Paper 1). In this section, we identify various mech- 
anisms and results from that paper, relevant to the current 
investigation, which have not been dealt with elsewhere. 

We use the same weight convergence assessment mech- 
anism as in Paper 1. It is important to note that if a M2M 
model reproduces the constraining observables, this does not 
necessarily mean that the particle weights have converged. 

We demonstrated in Paper 1 that using regularisation 
(/i 7^ 0) and luminosity density as a constraint made lit- 
tle difference to the model determined mass-to-light ratio. 
We have chosen not t o use regularisation in this paper. 
ICappellari et al.l (|2006h did in fact use integral space sec- 
ond derivative regularisation in their mass-to-light exercise. 
Similarly to regularisation, we choose not to use luminosity 
density as a constraint on the luminous matter distribution 
within our M2M models as a matter of course, but include 
it, where explicitly state d, for comparison purpo ses only. 

For consistency with ICappellari et al.l ([2006) , we quote 
no confidence intervals on our model determined mass-to- 
light ratios. 



4 MASS-TO-LIGHT DETERMINATION 

The process for determining the mass-to-light ratio for a 
galaxy is straightforward and widely used elsewhere. We run 
a series of M2M models varying a mass-to-light parameter 
(T) and look for a minimum in the resulting model x 2 values. 
The parameter value at the minimum, adjusted for the black 
hole mass, is taken as the 'true' mass-to-light ratio for the 
galaxy given all the modelling assumptions. 

T _ T m odel^ + MbH 

J- adjusted — ^~ V^V 

where L is the model total luminosity and Mbh is the black 
hole mass. 

The mass-to-light ratios we achieve are shown in Table 
[TJ and Figure [3] contains a plot of the M2M mass -to- lig ht val- 
ues against those achieved bv ICappellari et alJ (|2006l ) using 
Schwarzschild's method. For 96% of the galaxies, the relative 
difference between the mass-to-light values is < 20% and for 
83% of the galaxies, the difference is < 10%. 3 galaxies have 
differences > 15% - NGC 524 (28%), NGC 4486 (16%), and 
NGC5845 (17%). Table [2] contains a fuller breakdown. 

Performing a least squares straight line fit to the loga- 
rithmic mass-to-light data yields 

«r ^ ~a.ii±o.o4 / c\ 

1 M2M OC I gch (2,0) 

where Tm2M is the M2M mass-to-light ratio and Ts c h is the 
Schwarzschild equivalent. Introducing luminosity density to 
constrain the distribution of luminous matter does not al- 
ter this relationship. (For most galaxies, the models already 
well reproduce the density without the use of an explicit con- 
straint). Changing the particle initial conditions to scheme 
1 in section 13.51 (match the luminosity density spatially, ve- 
locities based on the Jeans' equations) does not significantly 
alter the relationship with 



0.8 
0.7 
0.6 
8 0.5 

E 0.4 

cn 
o 

0.3 

0.2 



0.3 0.4 0.5 0.6 

log (M/L) (Schwarzschild) 



Figure 3. Comparison of M2M vs Schwarzschild mass-to-light 
ratios with galaxies which differ by more than 10% marked in 
red. The solid blue line represents equality of the ratios, and the 
black dashed line a least squares fit to the data. 



Table 2. Comparison between M2M and SAURON mass-to-light 
ratios 



Difference 


Number of 
galaxies 




Galaxy 

NGC numbers 


< 5% 


9 


38% 


2974, 3377, 3608, 4150, 








4458, 4526, 4621, 5813, 








5846 


5% to 10% 


11 


46% 


821, 3156, 3379, 3414, 








4278, 4374, 4459, 4473, 








4550, 4552, 4660 


10% to 15% 


1 


4% 


7457 


> 15% 


3 


13% 


524, 4486, 5845 



oc To 



(26) 



Removing NGC 524 as an 'outlier' in the M2M results, we 
calculate the root mean square deviation of the two sets of 
mass-to-light ratios as 0.31. Assuming the errors in both the 
Schwarzschild and M2M methods to be similar, we arrive at 
an intrinsic error in the methods of 5.9% (calculated as the 
rms deviation divided by y/2 divided by the mean M2M 
mass-to-light r atio). This figure agree s well with the value 
(6%) quoted bv lCappellari et al.1 (|2006|) for their comparison 
of Jeans and Schwarzschild mode ls, and with the theo retical 
model values (« 5%) achieved in lLong &; Mad (|2010h . 

Despite differences in the two modelling methods, in 
general the methods are delivering, as one would hope, sim- 
ilar mass-to-light ratios for a galaxy. However, it is apparent 
from Figure [3] that, for the sample of galaxies analysed, ei- 
ther the M2M method is slightly over-estimating the mass- 
to-light ratios, or conversely, that Schwarzschild's method is 
under-estimating them. 

Having determined an estimate for the galaxies' mass- 
to-light ratios, we perform a further set of modelling runs 
at those ratios using 10 6 particles in order to investigate 
how well the observables are reproduced and to calculate 
the global anisotropy parameters (see section [5}. At the esti- 
mated mass-to-light ratios, weight convergence is good with 
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Table 3. x 2 per degree of freedom analysis - all galaxies 



Observable 


Minimum 

X 2 /dof 


Maximum 
X 2 /dof 


Mean 
X 2 /dof 


Median 
X 2 /dof 


sb 


0.05 


0.64 


0.30 


0.25 


hi 


0.05 


0.96 


0.35 


0.27 


h 2 


0.39 


1.22 


0.73 


0.75 


h 3 


0.19 


2.17 


0.58 


0.50 


/l4 


0.36 


1.20 


0.66 


0.61 


h 5 


0.000 


0.023 


0.004 


0.001 


he 


0.002 


0.021 


0.009 


0.008 


V m 


0.05 


0.97 


0.35 


0.27 


0~m 


0.33 


2.31 


0.78 


0.70 



The M2M constraining observables are surface brightness (sb) 
and hi to h6. The model mean line-of-sight velocity (%) and 
velocity dispersion (cr m ) are calculated as in equations 1161 and 
1171 The maximum values of hs and Vm come from the modelling 
of NGC 5845 and the maximum of cr m from NGC 3156. 

« 98% of particles having converged weights. The \ 2 P er 
degree of freedom values for the constraining observables 
(surface brightness and the Gauss-Hermite coefficients h\ to 
ho) are generally less than 1. The values for the calculated 
observables (mean line-of-sight velocity and velocity disper- 
sion) are similarly so. More detail is given in Table [3] In 
Figure 31 we show the observed and model line-of-sight ve- 
locity maps for a selection of 4 galaxies (NGC 2974, NGC 
3414, NGC 4550 and NGC 5813) covering the major galaxy 
types (elliptical and lenticular) , rotation (fast and slow) , ori- 
entation (inclined to the line-of-sight and edge-on) and core 
features. In particular, we note that the M2M method is 
able to reproduce kinematically distinct cores and counter 
rotating cores. For completeness, Figure [5] contains a com- 
plete set of observable maps (velocity, dispersion, /i3, and 
/14) for NGC 4660. Overall, reproduction is satisfactory. 



where T is the mass-to-light ratio of the galaxy, L is the 
model luminosity, akj is the mean luminosity weighted ve- 
locity dispersion in direction k in bin j, and Wj, the sum of 
the particles weights in bin j, is given by 



N 



Wi S(i £ j). 



(30) 



ICappellari et al. (2007) introduce two further parame- 
ters, P and 7. Using cylindrical polar coordinates (R,<f>,z), 



13 = 1 



and 



7 : 



(31) 



(32) 



The global anisotropy parameter 5 is then calculated as 

8=^. (33) 

As noted in ICappellari et aTl (|2007l ) . /3 describes the global 
shape of the velocity dispersion tensor in the (vr, v z ) plane, 
and 7 the shape in a plane orthogonal to v z . 

The three anisotropy parameters described so far, whilst 
still applicable to spherical galaxies, are more ap propriate 
to axisymmetric galaxies. ICappellari et all (J2007) describe 
a further parameter, /3 r , to measure the anisotropy of (near) 
spherical galaxies. 

Hoe + 



Pr = l- 



2 n r 



(34) 



where (r, (j)) are spherical coordinates. 

For all parameters, the radial and tangential velocity 
dispersion anisotropy regimes are shown in Figure 6. For 
the directions indicated within the definitions of the pa- 
rameters, a zero parameter value indicates isotropy and a 
positive value, a radial bias to the velocity dispersion. 



5 ANISOTROPY PARAMETERS 



5.2 M2M modelling 



5.1 Theory 

The global aniso tr opy parameter 5 is given in 
iBinnev fe Tremaind (|2008T ) §4.8 and ICappellari et all 
(2007) as 



1 



U x 



(27) 



where Uij is from the kinetic energy due to random stellar 
motion, z indicates the symmetry axis of an axisymmetric 
galaxy and x is some fixed direction orthogonal to it. Hkk is 
defined as 



n fc 



2 j3 

paf, a x 



(28) 



where p is the mass density and the velocity dispersion 
in direction k. The equivalent for M2M modelling purposes, 
calculated using the particle weights and binning particle 
data into J bins, is 



n fc 



2 



(29) 



For each galaxy, we perform a M2M modelling run using the 
mass-to- light ratio determined in section 2] We calculate the 
anisotropy parameters by binning the end of run particle 
velocity data on an (R,z) grid for f3 and 7, and for /3 r , 
we bin the data radially. Given we are using the particle 
data directly, the number of part i cles is increased to 1 x 
10 6 . Similarly to ICappellari et all |2007), we only include 
particles in the calculation which are currently within 25 
arcsec spherical radius of the galactic centre. 

Overall, we achieve reasonable (but not good) agree- 
ment with the SAURON values of the global anisotropy pa- 
rameter 5 with 14 of the 24 galaxies having M2M values 
equal to ^sauron ±0.05. The detailed results are captured 
in Table 2) Figure [6] sh ows /3, 7 and (3 r plotte d against 5 
and is the equivalent of ICappellari et all (|2007h Figure 2. 
Figures [7] and [8] compare the M2M parameter values against 
the SAURON values. Two main differences can be seen, the 
first being that in the M2M results the slow rotating galaxies 
are more clustered in parameter space. The second is that 
from Figure [6] the number of galaxies exhibiting tangential 
anisotropy has reduced and we examine this in more detail 
below. 
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Figure 4. Galaxy velocity maps for NGC 2974 (inclined, elliptical, fast rotator), NGC 3414 (edge-on, lenticular, slow rotator with a 
counter rotating core), NGC 4550 (inclined, lenticular, fast rotator with counter rotating disc) and NGC 5813 (edge-on, elliptical, slow 
rotator with a kinematically distinct core) showing the target mean line-of-sight velocity and the model produced versions. See Tabled 
for more data on the individual galaxies. The velocity units are km s — 1 . 
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Figure 5. The complete set of target and model observable maps for NGC 4660. 
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Figure 6. Anisotropy parameters (3, 7 and (3 r for the individual galaxi es plotted against S. Slo w rotating galaxies are plotted with red 
solid markers. The M2M plots are in the right hand column. Plots from Cappe lTari et all ([20071 ) (left hand column) have been included 
for comparison purposes. The black dashed lines separate the plot areas into radial (upper or upper left) and tangential (lower or lower 
right) anisotropy regimes. It is clear that fewer M2M models exhibit tangential velocity dispersion anisotropy and that the slow rotating 
galaxies are more tightly clustered in parameter space than their SAURON equivalents. 



We start by performing a simple count of the num- 
ber of positive and negative parameter values - the result 
is recorded in Table [5] For the fast rotating galaxies, any 
change of sign of a parameter between the SAURON and 
M2M values comes from one of 6 galaxies (NGC 2974, 4150, 



4278, 4473, 4621, 4660). Similarly, for slow rotating galax- 
ies, only 4 galaxies are involved (NGC 3414, 3608, 4458, 
4552). From examining the characteristics of the galaxies, 
there are no obvious groupings which might help explain 
the differences between the SAURON and M2M results. For 
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Table 4. Comparison between SAURON and M2M anisotropy parameters 



Cappellari et ah (2007) From M2M models 



Galaxy Inclination 
(cleg) 


Fast 
Rotator 




Pr 




p 




7 




6 




Pr 




p 




7 




s 


NGC 524 


19 


VPS 





06 





17 


-0 


04 





19 





00 





27 


-0 


21 





34 


NGC 821 


90 


VPS 


o 


16 


o 


21 


o 


04 


o 


20 


o 


25 


o 


25 


o 


13 


o 


20 


NGC 2974 


57 


VPS 


-0 


20 


o 


13 


-0 


30 


o 


24 


o 


11 


o 


13 


o 


03 


o 


12 


NGC 3156 


68 


VPS 


o 


17 


o 


39 


o 


19 


o 


33 


o 


24 


o 


19 


o 


33 


o 


03 


NGC 3377 


90 


VPS 

J 


o 


07 





28 


o 


08 


o 


25 





21 





22 





15 





16 


NGC 3379 


90 


VPS 


o 


11 


o 


06 


o 


06 


o 


03 


o 


33 


o 


15 


o 


26 


o 


02 


NGC 3414 


90 


no 


-0 


12 





06 


-0 


12 


o 


11 





25 





13 





15 





06 


NGC 3608 


90 


no 


o 


04 


o 


10 


-0 


06 


o 


13 


o 


23 


o 


17 


o 


10 


o 


13 


NGC 4150 


52 


VPS 


-0 


01 


o 


32 


-0 


12 


o 


36 


o 


15 


o 


06 


o 


17 


-0 


03 


NGC 4278 


90 


VPS 

J 


-0 


02 


o 


11 


-0 


17 


o 


18 


o 


33 


o 


29 


o 


21 


o 


20 


NGC 4374 


90 


no 


o 


11 


o 


10 


o 


05 


o 


08 


o 


30 


o 


17 


o 


21 


o 


07 


NGC 4458 


90 


no 


-0 


26 


-0 


01 


-0 


23 





09 





26 





13 





12 





07 


NGC 4459 


47 


yes 





10 





05 





11 





00 





23 





13 





16 





05 


NGC 4473 


73 


yes 


-0 


21 





18 


-0 


50 





34 





15 





25 


-0 


02 





26 


NGC 4486 


90 


no 





24 





11 





22 





00 





35 





15 





26 





03 


NGC 4526 


79 


yes 





11 





11 





09 





06 





21 





20 





16 





13 


NGC 4550 


84 


yes 


-0 


37 





43 


-0 


87 





60 


-0 


26 





27 


-0 


79 





48 


NGC 4552 


90 


no 


-0 


06 





01 


-0 


03 





02 





33 





16 





27 





03 


NGC 4621 


90 


yes 


-0 


04 





11 


-0 


17 





18 





17 





20 





03 





19 


NGC 4660 


70 


yes 





02 





27 


-0 


11 





30 





23 





27 





06 





25 


NGC 5813 


90 


no 





17 





18 





21 





08 





32 





18 





22 





08 


NGC 5845 


90 


yes 





24 





23 





18 





15 





36 





23 





31 





09 


NGC 5846 


90 


no 





17 





09 





17 





01 





36 





15 





25 





03 


NGC 7457 


64 


yes 





03 





38 





04 





37 





21 





21 





30 





07 



Comparison between the SAURON values for the anisotropy parameters p r , P, 7 and 5 taken from lCappellari et al.l (|2007l ) Table 2, and 
the same parameters calculated from M2M models. For the 8 parameter, 14 M2M values (including all the slow rotating galaxies) are 

within ±0.05 of the SAURON values. 



.•<s«V*' 



(3 (SAURON) 
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Y (SAURON) 
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Figure 8. Comparison between the M2M and SAURON determined anisotropy parameters P, 7 and fi r . Slow rotating galaxies are 
plotted with red solid markers. The black dashed lines indicate equality of the M2M and SAURON values. 



example, the 10 galaxies include both elliptical and lentic- 
ular galaxies and galaxies which are inclined to the line-of- 
sight or edge-on. All the galaxies have M2M mass- to- light 
ratios which differ from the SAURON values by < 7.5%, 
and 7 of the galaxies have M2M 5 anisotropy parameter 
values within ±0.05 of the SAURON values. It should be 
noted that for a given value of the global anisotropy param- 
eter 5, there is a linear relationship between f3 and 7 so a 
range of orbital models and velocity dispersion anisotropics 
is to be expected given the kinematic observations avail- 
able. Including luminosity density as a constraint does not 
alter th e tangential anis o tropy result. Note that the data 
used in I Cappellari et al.l (|2007h (see section 2) differ from 



the SAURON data release though it is not clear that this 
would explain the differences. 

In the absence of further constraints and perhaps more 
detailed investigations, we conclude that the differences in 
the anisotropy parameter values are due to differences in the 
particle / orbit initial conditions and the modelling methods 
used resulting in different orbital weightings. 



6 CONCLUSIONS 

We have undertaken the largest M2M exercise to date and 
re-analysed 24 elliptical and lenticular galaxies previously 
analysed with Schwarzschild's method. We have used the 
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Table 5. Anisotropy parameter analysis 



Number of Number of 

Positive Values Negative Values Galaxy 
Parameter Sauron M2M Sauron M2M NGC numbers 



Fast Rotators 

(3 r 10 15 6 1 2974, 4150, 4278, 4473, 4621 

P 16 16 

7 8 13 8 3 2974, 4150, 4278, 4621, 4660 

S 16 15 1 4150 
Slow Rotators 

(3 r 5 8 3 3414, 4458, 4552 

(3 7 8 1 4458 

7 4 8 4 3414, 3608, 4458, 4552 

S 8 8 



The table shows the numbers of galaxies with positive and negative velocity dispersion anisotropy parameter values as a means of 
assessing any general change of radial or tangential anisotropy regime for the sample of galaxies. The 'galaxy' column indicates those 

galaxies where the regime differs between the SAURON and M2M models. 







o 




, ^ o 




Q • , o 
^ *• ° 

o 





"°-0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 



5 (SAURON) 

Figure 7. Comparison between the M2M and SAURON de- 
termined values of the global anisotropy parameter 5. Slow ro- 
tating galaxies are plotted with red solid markers. The black 
dotted lines indicate ± 0.05 of the SAURON values recorded in 
ICappellari et al.l (|2007T ). 

M2M method as far as possible as a 'black box' - we have de- 
veloped no computer code specific to any one galaxy. Where 
there are modelling parameters to be set or tuned, we have 
adopted the same strategy and process for all galaxies. 

Our M2M implementation has been adapted to use ob- 
servable data available as a Voronoi tessellation and to han- 
dle gravitational potentials derived from a mult i- Gaussian 
expansion and deprojection of a galaxy's surface brightness. 
We have identified a computer performance issue with our 
M2M implementation which may affect other users of the 
method depending on their implementation and network 
configuration. For the future, an improved process, prefer- 
ably computerised, for setting the global rotation of the sys- 
tem of particles is required. 

We achieve reasonable agreement (14 out of 24 galaxies) 
with the SAURON values of the global anisotropy parame- 
ter S but our overall assessment is that further (theoretical) 
investigations of the impact of orbit / particle initial condi- 
tions and the resultant orbit weights are required before the 



differences between the SAURON and M2M methods can be 
fully understood. In the M2M case, it may prove to better 
to calculate an exponentially smoothed version of 5 rather 
than relying on the end of modelling run particle data. 

We have demonstrated that, despite differences in the 
M2M and Schwarzschild modelling methods, in general the 
methods are delivering similar mass-to-light ratios for a 
galaxy. Whether the slight over estimation (M2M) or un- 
der estimation (Schwarzschild) is a real effect or not will 
only be resolved by using a different sample of galaxies. 
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